## A file assess the relationships between
## perceptions and variables representing alternative explanation

options(digits = 4, width = 132)

library(here)
library(xtable)
library(lmerTest)
library(tidyverse)
library(sf)

load(here("Analysis", "analysis_anyDA_new.rda"), verbose = TRUE)
load(here("Data", "wrkdatOwnMap_new.rda"), verbose = TRUE)
load(here("Design", "matches_anyDA_new.rda"), verbose = TRUE)

table(is.na(wrkdatOwnMap_new$vm.community.norm2))

wdat0 <- inner_join(wrkdatOwnMap_new %>% select(-geometry), matches_anyDA_new)

## Add information about which unit is ranked higher versus lower within pair
wdat1 <- wdat0 %>%
  group_by(pair) %>%
  mutate(
    perc_more = rank(vm.community.subj) - 1,
    social_capital01_rank = rank(social.capital) - 1,
    community_resp01_rank = rank(community.resp01) - 1,
    s_racism_rank = rank(s.racism) - 1,
    norm2_rank = rank(vm.community.norm2) - 1
  ) %>%
  ungroup() %>%
  select(
    vcid, dauid, pair, vm.community.subj, vm.community.norm2, perc_more, norm2_rank,
    social.capital01, community.resp01,
    da_prop_vm_20pct_06, s_racism_rank, s.racism,
    racial.conflict.work.ethnicity, minorities.community.ethnic.friends.composition,
    minorities.community.group.feeling.thermometers.other.asian,
    minorities.community.group.feeling.thermometers.east.indian, minorities.community.group.feeling.thermometers.chinese,
    minorities.community.group.feeling.thermometers.black, minorities.community.group.feeling.thermometers.latin,
    minorities.community.group.feeling.thermometers.white, educnew,
    social_capital01_rank, community_resp01_rank
  ) %>%
  droplevels()

## We want to look at whether s.racism -> perceptions AND s.racism -> the
## different outcomes. It might be clearer to ask whether the person higher in
## s.racism is also the person higher in.norm2ective perceptions.

stopifnot(all(wdat1$s_racism_rank >= 0 & wdat1$s_racism_rank <= 1))

## Create ranked versions of both work contact and friend contact variables

## Are the people at your work (or school) mostly white, mostly racial or
## ethnic minorities, about half and half, or some other mixture?

## Now thinking more generally, are your friends mostly white, mostly racial
## and ethnic minorities, about half and half, or some other mixture of people?


wdat1 <- wdat1 %>%
  mutate(
    workcontact01 = case_when(
      racial.conflict.work.ethnicity == "white" ~ 0,
      racial.conflict.work.ethnicity == "other" ~ 1,
      racial.conflict.work.ethnicity == "ethnic" ~ 1,
      racial.conflict.work.ethnicity == "half" ~ 1
    ),
    friendcontact01 = case_when(
      minorities.community.ethnic.friends.composition == "white" ~ 0,
      minorities.community.ethnic.friends.composition == "other" ~ 1,
      minorities.community.ethnic.friends.composition == "ethnic" ~ 1,
      minorities.community.ethnic.friends.composition == "half" ~ 1
    ),
    relativeFT = (minorities.community.group.feeling.thermometers.other.asian
    + minorities.community.group.feeling.thermometers.east.indian
      + minorities.community.group.feeling.thermometers.chinese
      + minorities.community.group.feeling.thermometers.black
      + minorities.community.group.feeling.thermometers.latin) / 5
      - minorities.community.group.feeling.thermometers.white,
  ) %>%
  group_by(pair) %>%
  mutate(
    workcontact01_ranked = rank(workcontact01) - 1,
    friendcontact01_ranked = rank(friendcontact01) - 1,
    relativeFT_ranked = rank(relativeFT) - 1
  ) %>%
  ungroup()


### Checking to see if the code acheived what we want
## The higher value is 1 and the lower is 0
wdat1 %>%
  filter(pair %in% unique(wdat1$pair)[1:3]) %>%
  group_by(pair) %>%
  select(s.racism, s_racism_rank) %>%
  arrange(pair, s.racism)

with(wdat1, table(s_racism_rank, s.racism, exclude = c()))
with(wdat1, table(perc_more, vm.community.norm2, exclude = c()))
with(wdat1, table(perc_more, vm.community.subj, exclude = c()))
with(wdat1, table(s_racism_rank, perc_more, exclude = c()))

## What direction is s.racism coded? This is just a double check  given the
## lack of a formal code book.

wdat1$college_plus <- with(wdat1, ifelse(educnew %in% c("bachelor's degree", "professional degree or doctorate", "master's degree"), 1, 0))
with(wdat1, table(college_plus, educnew, exclude = c()))

with(wdat1, table(college_plus, s.racism))
lm(s.racism ~ college_plus, data = wdat1)

## So higher values on s.racisms are more racist because of what we know about
## education and responses to the symbolic racism scale

##################### The analyses ############
### Do the analyses and interpret each result. What does it tell us about how
### s.racism matters for these variables?

####
# vm.community.norm2
# social.capital01
# community.resp01

## Within-pair comparison on perception
wdat1 %>%
  group_by(pair) %>%
  mutate(
    perc_diff2 = abs(diff(vm.community.norm2)),
    perc_diff = abs(diff(vm.community.subj))
  ) %>%
  ungroup() %>%
  distinct(pair, .keep_all = TRUE) %>%
  summarise(
    Mean = mean(perc_diff, na.rm = TRUE),
    SD = sd(perc_diff, na.rm = TRUE),
    "%25 Q" = quantile(perc_diff, .25, na.rm = TRUE),
    "%50 Q" = quantile(perc_diff, .5, na.rm = TRUE),
    "%75 Q" = quantile(perc_diff, .75, na.rm = TRUE),
    Mean2 = mean(perc_diff2, na.rm = TRUE),
    SD2 = sd(perc_diff2, na.rm = TRUE),
    "%25 Q2" = quantile(perc_diff2, .25, na.rm = TRUE),
    "%50 Q2" = quantile(perc_diff2, .5, na.rm = TRUE),
    "%75 Q2" = quantile(perc_diff2, .75, na.rm = TRUE)
  )
# A tibble: 1 × 10
#    Mean    SD `%25 Q` `%50 Q` `%75 Q` Mean2   SD2 `%25 Q2` `%50 Q2` `%75 Q2`
#   <dbl> <dbl>   <dbl>   <dbl>   <dbl> <dbl> <dbl>    <dbl>    <dbl>    <dbl>
# 1 0.248 0.263  0.0800    0.17    0.34 0.228 0.201     0.07     0.17     0.33

## A check on the code below:
pair_summaries <- wdat1 %>%
  group_by(pair) %>%
  dplyr::select(workcontact01, friendcontact01) %>%
  summarize(across(where(is.numeric), ~ abs(diff(.x))))

summary_function <- function(x) {
  x <- x[!is.na(x)]
  c(mean = mean(x), sd = sd(x), quantile(x, c(.1, .25, .5, .75, .9)))
}

pair_summaries %>%
  reframe(across(where(is.numeric), summary_function)) %>%
  mutate(stat = c("mean", "sd", "q10", "q25", "q50", "q75", "q90"))


## Within-pair comparison on work contact
wdat1 %>%
  group_by(pair) %>%
  mutate(perc_diff = abs(diff(workcontact01))) %>%
  ungroup() %>%
  distinct(pair, .keep_all = TRUE) %>%
  summarise(
    Mean = mean(perc_diff, na.rm = TRUE),
    SD = sd(perc_diff, na.rm = TRUE),
    "%25 Q" = quantile(perc_diff, .25, na.rm = TRUE),
    "%50 Q" = quantile(perc_diff, .5, na.rm = TRUE),
    "%75 Q" = quantile(perc_diff, .75, na.rm = TRUE)
  )

## A tibble: 1 × 5
#   Mean    SD `%25 Q` `%50 Q` `%75 Q`
#  <dbl> <dbl>   <dbl>   <dbl>   <dbl>
# 1 0.433 0.496       0       0       1

## Within-pair comparison on friend contact
wdat1 %>%
  group_by(pair) %>%
  mutate(perc_diff = abs(diff(friendcontact01))) %>%
  ungroup() %>%
  distinct(pair, .keep_all = TRUE) %>%
  summarise(
    Mean = mean(perc_diff, na.rm = TRUE),
    SD = sd(perc_diff, na.rm = TRUE),
    "%25 Q" = quantile(perc_diff, .25, na.rm = TRUE),
    "%50 Q" = quantile(perc_diff, .5, na.rm = TRUE),
    "%75 Q" = quantile(perc_diff, .75, na.rm = TRUE)
  )

# %# A tibble: 1 × 5
# %   Mean    SD `%25 Q` `%50 Q` `%75 Q`
# %  <dbl> <dbl>   <dbl>   <dbl>   <dbl>
# %1 0.351 0.477       0       0       1

## work contact
wdat1 %>%
  filter(workcontact01_ranked != 0.5) %>%
  group_by(pair) %>%
  filter(n() > 1) %>%
  arrange(workcontact01_ranked) %>%
  mutate(diff = vm.community.norm2 - lag(vm.community.norm2)) %>%
  ungroup() %>%
  summarise(mean = mean(diff, na.rm = TRUE), sd = sd(diff, na.rm = TRUE))

## A tibble: 1 × 2
#    mean    sd
#   <dbl> <dbl>
# 1 0.0115 0.300


## friend contact
wdat1 %>%
  filter(friendcontact01_ranked != 0.5) %>%
  group_by(pair) %>%
  filter(n() > 1) %>%
  arrange(friendcontact01_ranked) %>%
  mutate(diff = vm.community.norm2 - lag(vm.community.norm2)) %>%
  ungroup() %>%
  summarise(mean = mean(diff, na.rm = TRUE), sd = sd(diff, na.rm = TRUE))
# %# A tibble: 1 × 2
# %    mean    sd
# %   <dbl> <dbl>
# %1 0.0626 0.316

## relative FT
wdat1 %>%
  filter(relativeFT_ranked != 0.5) %>%
  group_by(pair) %>%
  filter(n() > 1) %>%
  arrange(relativeFT_ranked) %>%
  mutate(diff = vm.community.norm2 - lag(vm.community.norm2)) %>%
  ungroup() %>%
  summarise(mean = mean(diff, na.rm = TRUE), sd = sd(diff, na.rm = TRUE))

## A tibble: 1 × 2
#     mean    sd
#    <dbl> <dbl>
# 1 -0.0114 0.305

## Create mean DA VM within each pair.
wdat1 <- wdat1 %>%
  group_by(pair) %>%
  mutate(da_prop_vm_20pct_06bi = mean(da_prop_vm_20pct_06)) %>%
  ungroup()
wdat1$da_prop_vm_20pct_06bi2 <- ave(wdat1$da_prop_vm_20pct_06, wdat1$pair)
all.equal(wdat1$da_prop_vm_20pct_06bi, wdat1$da_prop_vm_20pct_06bi2)

sr_perc <- lmer(vm.community.norm2 ~ s.racism + da_prop_vm_20pct_06 + (1 | pair) + (1 | dauid), data = wdat1)

sr_perc_rank <- lmer(norm2_rank ~ s_racism_rank + da_prop_vm_20pct_06 + (1 | dauid), data = filter(wdat1, norm2_rank != .5 & s_racism_rank != .5))

allfits_1 <- allFit(sr_perc_rank)
summary(allfits_1)$which.OK
## Result: All True
summary(sr_perc)$coef
#                     Estimate Std. Error   df t value   Pr(>|t|)
# (Intercept)          0.25376   0.009781 3600  25.943 3.449e-136
# s.racism            -0.06088   0.016687 3691  -3.648  2.676e-04
# da_prop_vm_20pct_06  0.99867   0.034782 1887  28.712 9.836e-151
summary(sr_perc_rank)$coef
#                         Estimate Std. Error   df     t value   Pr(>|t|)
# (Intercept)          0.525821304    0.01335 3405 39.38811385 6.208e-280
# s_racism_rank       -0.051643192    0.01711 3405 -3.01752967  2.567e-03
# da_prop_vm_20pct_06  0.000003808    0.07341 3405  0.00005187  1.000e+00

sr_soc_cap <- lmer(social.capital01 ~ s.racism + da_prop_vm_20pct_06 + (1 | pair) + (1 | dauid), data = wdat1)
sr_soc_cap_rank <- lmer(social_capital01_rank ~ s_racism_rank + da_prop_vm_20pct_06 + (1 | dauid),
  data = filter(wdat1, s_racism_rank != .5 & social_capital01_rank != .5)
)
allfits_2 <- allFit(sr_soc_cap_rank)
summary(allfits_2)$which.OK
## Result: All True
fixef(sr_soc_cap)
#         (Intercept)            s.racism da_prop_vm_20pct_06
#             0.72294            -0.01257            -0.15726
fixef(sr_soc_cap_rank)
#         (Intercept)       s_racism_rank da_prop_vm_20pct_06
#          0.51536677         -0.03073624          0.00001759


comm_resp <- lmer(community.resp01 ~ s.racism + da_prop_vm_20pct_06 + (1 | pair) + (1 | dauid), data = wdat1)
comm_resp_rank <- lmer(community_resp01_rank ~ s_racism_rank + da_prop_vm_20pct_06 + (1 | dauid),
  data = filter(wdat1, s_racism_rank != .5 & community_resp01_rank != .5)
)
allfits_3 <- allFit(comm_resp_rank)
summary(allfits_3)$which.OK
fixef(comm_resp)
#         (Intercept)            s.racism da_prop_vm_20pct_06
#             0.78117            -0.04868            -0.21575
fixef(comm_resp_rank)
#         (Intercept)       s_racism_rank da_prop_vm_20pct_06
#           0.5496670          -0.0993377           0.0000244

save(sr_perc, sr_perc_rank, sr_soc_cap, sr_soc_cap_rank, comm_resp,
  comm_resp_rank,
  file = here("Analysis", "alt_explanations_analysis.rda")
)
